Importar librerías¶

  • Pandas para el uso de datos
  • NumPy para operaciones matemáticas
  • Scipy para ajuste de curvas
  • Matplotlib.pyplot para graficar
  • Sklearn.metrics para calcular métricas
  • SQLAlchemy para conexión a base de datos
In [21]:
import pandas as pd
import numpy as np
from scipy.optimize import curve_fit
from scipy.optimize import minimize
from scipy.optimize import least_squares
import matplotlib.pyplot as plt
from sklearn.metrics import mean_squared_error, r2_score, mean_absolute_error
import scipy 
from sqlalchemy import create_engine, inspect
from sqlalchemy.orm import sessionmaker
import datetime

Conexión a base de datos¶

Para la base de datos se utiliza el servicio de https://neon.tech/, el cual se utiliza por integrar PostGIS sin costo y ofrecer una capacidad suficiente para pruebas.

In [16]:
db_url = "postgresql://test_owner:lpPrcHty8i9V@ep-royal-glitter-a50h4o6a.us-east-2.aws.neon.tech/test?sslmode=require"
engine = create_engine(db_url)

Lectura de Base de datos¶

Se guardó toda la información de la Especie con identificador 1 para hacer pruebas.

In [5]:
query = "SELECT * FROM data_especie1;"
df = pd.read_sql_query(query, engine)
df
Out[5]:
PARCELA # ROTACIÓN # ZONA AÑO PLANTACIÓN EDAD MEDICIÓN # FECHA MEDICIÓN ESPECIE PROCEDENCIA GENÉTICA ARBOL # ... DAP MEDIO [cm] AREA BASAL [m2/ha] H TOTAL [m] H DOMINANTE [m] H PODA [m] N/ha PODADOS [#] VOLUMEN EN PIÉ d.u. = 5 cm. [m3/ha scc]. IMA s/raleo [m3/ha/año]. Observación Clasificacion Densidad
0 1 1 1004 2012 6.043836 2 2018-06-20 00:00:00 1 Livingston HSCA T416 1 ... 17.229788 14.112918 9.768085 10.925 3.555319 None 51.132477 8.460269 None 3
1 1 1 1004 2012 6.043836 2 2018-06-20 00:00:00 1 Livingston HSCA T416 2 ... 17.229788 14.112918 9.768085 10.925 3.555319 None 51.132477 8.460269 None 3
2 1 1 1004 2012 6.043836 2 2018-06-20 00:00:00 1 Livingston HSCA T416 3 ... 17.229788 14.112918 9.768085 10.925 3.555319 None 51.132477 8.460269 None 3
3 1 1 1004 2012 6.043836 2 2018-06-20 00:00:00 1 Livingston HSCA T416 4 ... 17.229788 14.112918 9.768085 10.925 3.555319 None 51.132477 8.460269 None 3
4 1 1 1004 2012 6.043836 2 2018-06-20 00:00:00 1 Livingston HSCA T416 5 ... 17.229788 14.112918 9.768085 10.925 3.555319 None 51.132477 8.460269 None 3
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
383240 106 1 1004 2012 10.909589 7 2023-07-10 00:00:00 1 HSCB R 35 ... 28.300000 25.390000 17.600000 18.000 4.900000 None 173.100000 15.866775 None 2
383241 106 1 1004 2012 10.909589 7 2023-07-10 00:00:00 1 HSCB R 36 ... 28.300000 25.390000 17.600000 18.000 4.900000 None 173.100000 15.866775 None 2
383242 106 1 1004 2012 10.909589 7 2023-07-10 00:00:00 1 HSCB R 37 ... 28.300000 25.390000 17.600000 18.000 4.900000 None 173.100000 15.866775 None 2
383243 106 1 1004 2012 10.909589 7 2023-07-10 00:00:00 1 HSCB R 38 ... 28.300000 25.390000 17.600000 18.000 4.900000 None 173.100000 15.866775 None 2
383244 106 1 1004 2012 10.909589 7 2023-07-10 00:00:00 1 HSCB R 39 ... 28.300000 25.390000 17.600000 18.000 4.900000 None 173.100000 15.866775 None 2

383245 rows × 42 columns

Definición de modelos¶

Para el crecimiento se definieron 12 modelos de crecimiento. El ejercicio consiste en encontrar, a partir de los datos, los mejores valores constantes para cada modelo.

In [11]:
#Sweda
#N(t) = N0 * e^(r(t-t0))
def SWEDA(x,N0,r,t0):
    return N0 * np.exp(r*(x - t0))
    #return x * N0
#La fórmula para el modelo de crecimiento de Gompertz es la siguiente:
#N(t) = N0 * e^(-a * e^(-bt))
def GOMPERTZ(x,N0,a,b):
    return N0 * np.exp(-a * np.exp(-b * x))

#La función logística se puede escribir de la siguiente manera:#
#N(t) = K / (1 + ae^(-rt))
def LOGISTICA(x,k,a,r):
    return k / (1 + a * np.exp(-r * x))

def SIGMOID(x, L, k ,x0):
    """
    Modelo sigmoidal (función logística).

    Parámetros:
    - x: Array de valores independientes.
    - L: Límite superior de la curva.
    - x0: Punto medio de la curva en el eje x.
    - k: Pendiente de la curva en el punto medio.

    Retorna:
    - Array con los valores de la función sigmoidal para cada valor de x.
    """
    return L / (1 + np.exp(-k * (x - x0)))
def GOMPERTZ(x,N0,a,b):
    return N0 * np.exp(-a * np.exp(-b * x))

def GOMPERTZ2(x,N0,p,b):
    return N0 * np.exp(- np.exp(p -b * x))

def relacionPolimorfica(x,a,b,c):
    return a * ( 1 - np.exp( -b * x)) ** c

def MITSCHERLICH(x,a,L,b):
    return a * ( 1 - L* np.exp( -b * x))

#
def ChapmanRichards(x,a,b,c):
    return a * (1 -np.exp( -b * x)) ** c + np.exp(1)
                
def HossfeldI(x,a,b,c):
    return ((x ** 2) / (a + b*x) ** 2) + np.exp(1)
                
def Schumacher(x,a,b,c):
    return a * np.exp(-b * ( 1/ x)) + np.exp(1)
                
def Weibull(x,a,b,c):
    return a * ( 1 - np.exp( -b * x ** c)) + np.exp(1)

Funciones para modelar¶

En esta sección se definen un conjunto de funciones para encontrar el mejor modelo:

  • evaluar_modelo: Genera diferentes métricas para los valores encontrados.
  • getMinimosCuadrados: A través de la técnica de mínimos cuadrados de manera iterativa se encuentran los mejores valores para ajustar los modelos de la mejor manera.
  • getGrafico2: Con esta función se grafican los datos y la curva ajustada.
  • MODELOS2: Se evalúan todos los modelos definidos anteriormente.
In [31]:
def evaluar_modelo(funcion, X, Y,popt):  
    predicciones = funcion(X, *popt)
    # Calcular métricas básicas
    mse = mean_squared_error(Y, predicciones)
    r2 = r2_score(Y, predicciones)
    mae = mean_absolute_error(Y, predicciones)
    mape = np.mean(np.abs((Y - predicciones) / Y)) * 100
    std_residuos = np.std(Y - predicciones)
    
    # Otras métricas adicionales
    suma_residuos_abs = np.sum(np.abs(Y - predicciones))
    max_error_abs = np.max(np.abs(Y - predicciones))
    mean_error_abs = np.mean(np.abs(Y - predicciones))
    
    # Métricas adicionales
    suma_residuos_cuad = np.sum((Y - predicciones)**2)
    mediana_residuos = np.median(Y - predicciones)
    percentil_90 = np.percentile(Y - predicciones, 90)
    varianza_residuos = np.var(Y - predicciones)
    suma_cuad_residual_abs = np.sum(np.abs((Y - predicciones)) ** 2)
    min_error_abs = np.min(np.abs(Y - predicciones))
    skew_residuos = scipy.stats.skew(Y - predicciones)
    kurtosis_residuos = scipy.stats.kurtosis(Y - predicciones)
    coef_variacion_residuos = np.std(Y - predicciones) / np.mean(Y - predicciones)
    error_cuadratico_medio = np.sqrt(np.mean((Y - predicciones) ** 2))
    
    metricas = {
        "MSE": mse,
        "R2": r2,
        "MAE": mae,
        "MAPE": mape,
        "STD Residuos": std_residuos,
        "Suma Residuos Abs": suma_residuos_abs,
        "Max Error Abs": max_error_abs,
        "Mean Error Abs": mean_error_abs,
        "Suma Residuos Cuadráticos": suma_residuos_cuad,
        "Mediana Residuos": mediana_residuos,
        "Percentil 90": percentil_90,
        "Varianza Residuos": varianza_residuos,
        "Suma Cuadrada de los Residuos Absolutos": suma_cuad_residual_abs,
        "Min Error Abs": min_error_abs,
        "Skewness de los Residuos": skew_residuos,
        "Kurtosis de los Residuos": kurtosis_residuos,
        "Coeficiente de Variación de los Residuos": coef_variacion_residuos,
        "Error Cuadrático Medio": error_cuadratico_medio
    }

    return metricas


def getMinimosCuadrados(funcion,X,Y):
    def residuals(params, x, y):
        return y - funcion(x, *params)
    initial_guess = [1,1,1]
    result = least_squares(residuals, initial_guess, args=(X, Y))
    N0, t0, r = result.x
    return result.x

# Establecer la visualización en línea y mejorar la calidad de la imagen
%matplotlib inline
%config InlineBackend.figure_format = 'retina'

def getGrafico2(funcion, X, Y):
    # Generar puntos para el ajuste
    t_fit = np.linspace(min(X), max(X), 100) 
    
    # Obtener los parámetros del ajuste
    a,b,c = getMinimosCuadrados(funcion,X,Y)
    
    print(funcion.__name__,a,b,c)
    # Evaluar la función ajustada en los puntos de ajuste
    metricas = evaluar_modelo(funcion, X, Y,[a,b,c])
    N_fit = funcion(t_fit,a,b,c)
    
    # Plotear los datos y el ajuste
    plt.figure(figsize=(15, 9))  # Tamaño de la figura
    plt.scatter(X, Y, color='blue', label='Datos')  # Plotear los datos
    plt.plot(t_fit, N_fit, 'r-', label=f'Ajuste {funcion.__name__}')  # Plotear el ajuste
    plt.xlabel('Tiempo [Años]')
    plt.ylabel('DAP [mm]')
    plt.legend()
    plt.grid(True)  # Agregar una rejilla
    plt.title('Gráfico de Ajuste')  # Título del gráfico
    plt.tight_layout()  # Ajustar el diseño de la figura para evitar superposiciones
    plt.show()
    
    # Obtener métricas del modelo
    #metricas = evaluar_modelo(funcion, X, Y, params)
    metricas["Funcion"] = funcion.__name__
    metricas["b0"] = a
    metricas["b1"] = b
    metricas["b2"] = c
    return metricas

def MODELOS2(X,Y):
    d1 = getGrafico2(SWEDA,X,Y)
    d2 = getGrafico2(GOMPERTZ,X,Y)
    d3 = getGrafico2(LOGISTICA,X,Y)
    d4 = getGrafico2(SIGMOID,X,Y)
    d5 = getGrafico2(GOMPERTZ2,X,Y)
    d6 = getGrafico2(relacionPolimorfica,X,Y)
    d7 = getGrafico2(MITSCHERLICH,X,Y)
    d8 = getGrafico2(ChapmanRichards,X,Y)
    d9 = getGrafico2(HossfeldI,X,Y)
    d10 =getGrafico2(Schumacher,X,Y)
    d11 =getGrafico2(Weibull,X,Y)
    return  [d1,d2,d3,d4,d5,d6,d7,d8,d9,d10,d11]
In [29]:
dfClean = df[['Clasificacion Densidad','EDAD','DAP [mm]']].dropna()

Ejecución de todos los modelos¶

Esta es la parte del proceso más larga y que requiere mayor capacidad de cómputo.

In [32]:
acumulador = []
for i in dfClean['Clasificacion Densidad'].sort_values().unique():
    print(i)
    aux2 = dfClean[dfClean['Clasificacion Densidad'] == i] 
    auxX = np.array(aux2['EDAD'])
    auxY = np.array(aux2['DAP [mm]'])
    
    data = pd.DataFrame(MODELOS2(auxX,auxY))
    data["Clasificacion Densidad"] = i
    acumulador.append(data.copy())
    break
errores = pd.concat(acumulador)
errores["Fecha_Analisis"] = datetime.datetime.now()
1
SWEDA 21.04119438647293 0.046419999139448886 7.10587799699723
GOMPERTZ 38.46327179557505 1.725621174054583 0.13785611393744854
LOGISTICA 37.032921601250976 3.2109256850486454 0.18705747048150667
SIGMOID 37.03315033104063 0.18705474512879688 6.236430348041789
GOMPERTZ2 38.46292029061756 0.5455947007342655 0.13785926728915074
C:\Users\limc_\AppData\Local\Temp\ipykernel_25800\2818810145.py:37: RuntimeWarning: invalid value encountered in power
  return a * ( 1 - np.exp( -b * x)) ** c
relacionPolimorfica 44.160254594609384 0.0685756872552295 0.8134078716935539
MITSCHERLICH 41.366632283607444 0.9515908033103654 0.08787302045324784
C:\Users\limc_\AppData\Local\Temp\ipykernel_25800\2818810145.py:44: RuntimeWarning: invalid value encountered in power
  return a * (1 -np.exp( -b * x)) ** c + np.exp(1)
ChapmanRichards 40.27005479755286 0.07826878248230645 0.9718266320152252
HossfeldI -0.6570581261937979 -0.1451961919079237 1.0
Schumacher 41.20107176680019 5.904595978915142 1.0
Weibull 41.09688061447227 0.08410061768487219 0.9664273943593613
In [35]:
errores
Out[35]:
MSE R2 MAE MAPE STD Residuos Suma Residuos Abs Max Error Abs Mean Error Abs Suma Residuos Cuadráticos Mediana Residuos ... Skewness de los Residuos Kurtosis de los Residuos Coeficiente de Variación de los Residuos Error Cuadrático Medio Funcion b0 b1 b2 Clasificacion Densidad Fecha_Analisis
0 27.558134 0.468702 4.029327 17.917199 5.249462 12434.504221 25.503244 4.029327 85044.402026 -0.300545 ... 0.368386 1.467820 -1.467915e+02 5.249584 SWEDA 21.041194 0.046420 7.105878 1 2024-05-02 16:56:46.900983
1 24.920378 0.519556 3.803969 16.033188 4.992028 11739.048860 24.255991 3.803969 76904.285665 -0.309014 ... 0.386523 1.732474 -8.306736e+02 4.992031 GOMPERTZ 38.463272 1.725621 0.137856 1 2024-05-02 16:56:46.900983
2 25.099417 0.516104 3.824390 16.227436 5.009920 11802.066760 24.268368 3.824390 77456.799770 -0.294159 ... 0.381045 1.697224 -4.531408e+02 5.009932 LOGISTICA 37.032922 3.210926 0.187057 1 2024-05-02 16:56:46.900983
3 25.099417 0.516104 3.824391 16.227448 5.009920 11802.069612 24.268388 3.824391 77456.799764 -0.294147 ... 0.381043 1.697223 -4.531715e+02 5.009932 SIGMOID 37.033150 0.187055 6.236430 1 2024-05-02 16:56:46.900983
4 24.920378 0.519556 3.803968 16.033172 4.992028 11739.045338 24.255969 3.803968 76904.285673 -0.308969 ... 0.386525 1.732475 -8.306188e+02 4.992031 GOMPERTZ2 38.462920 0.545595 0.137859 1 2024-05-02 16:56:46.900983
5 24.692280 0.523954 3.777398 15.808899 4.969129 11657.049170 24.326984 3.777398 76200.377269 -0.262152 ... 0.398153 1.787827 -8.413321e+02 4.969133 relacionPolimorfica 44.160255 0.068576 0.813408 1 2024-05-02 16:56:46.900983
6 24.727972 0.523266 3.779970 15.795186 4.972723 11664.988869 24.263190 3.779970 76310.521685 -0.236278 ... 0.391794 1.774221 7.188269e+06 4.972723 MITSCHERLICH 41.366632 0.951591 0.087873 1 2024-05-02 16:56:46.900983
7 24.730752 0.523212 3.782651 15.864034 4.972995 11673.260247 24.324262 3.782651 76319.100526 -0.287542 ... 0.397526 1.779102 -5.878602e+02 4.973002 ChapmanRichards 40.270055 0.078269 0.971827 1 2024-05-02 16:56:46.900983
8 24.684556 0.524103 3.771243 15.669904 4.968355 11638.056309 24.195639 3.771243 76176.539028 -0.198491 ... 0.401314 1.786289 2.991313e+03 4.968355 HossfeldI -0.657058 -0.145196 1.000000 1 2024-05-02 16:56:46.900983
9 24.873829 0.520454 3.782467 15.595337 4.987348 11672.694631 24.035584 3.782467 76760.637141 -0.173106 ... 0.413218 1.758948 3.638694e+02 4.987367 Schumacher 41.201072 5.904596 1.000000 1 2024-05-02 16:56:46.900983
10 24.728393 0.523257 3.782949 15.876494 4.972754 11674.181161 24.340027 3.782949 76311.819390 -0.284092 ... 0.399331 1.781259 -4.749267e+02 4.972765 Weibull 41.096881 0.084101 0.966427 1 2024-05-02 16:56:46.900983

11 rows × 24 columns

Guardar la base de datos de errores y constantes¶

In [26]:
errores.to_sql("data_error", engine, if_exists='append', index=False)
Out[26]:
55
In [ ]: